function [correlation_results, coherence_params, fig_correlation] = time_frequency_correlation(carrier_freq, fs, mobile_speed, delay_spread)
% 时频相关性分析
% 输入参数：
%   carrier_freq - 载波频率 (Hz)
%   fs - 采样频率 (Hz)
%   mobile_speed - 移动速度 (m/s)
%   delay_spread - 时延扩展 (秒)
% 输出参数：
%   correlation_results - 相关性结果
%   coherence_params - 相干参数
%   fig_correlation - 图形句柄

% 添加路径
addpath('../Common');
colors = color_definitions();

%% 并行计算设置
[use_parallel, pool_info] = parallel_manager();

% 参数设置
c = 3e8; % 光速

% 计算基本参数
fd = mobile_speed * carrier_freq / c; % 最大多普勒频移
fd_normalized = fd / fs; % 归一化多普勒频移

fprintf('最大多普勒频移: %.1f Hz\n', fd);
fprintf('归一化多普勒频移: %.4f\n', fd_normalized);

%% 时间相关性分析
tau = 0:0.001:0.1; % 时间间隔 (秒)

if use_parallel
    % 使用并行计算计算时间相关性
    fprintf('使用并行计算计算时间相关性...\n');
    parfor i = 1:length(tau)
        % Jakes模型 (经典多普勒谱)
        jakes_corr(i) = besselj(0, 2*pi*fd*tau(i));
        % 均匀多普勒谱
        uniform_corr(i) = sinc(2*fd*tau(i));
        % 高斯多普勒谱
        gaussian_corr(i) = exp(-(pi*fd*tau(i)).^2);
    end
else
    % Jakes模型 (经典多普勒谱)
    jakes_corr = besselj(0, 2*pi*fd*tau);
    % 均匀多普勒谱
    uniform_corr = sinc(2*fd*tau);
    % 高斯多普勒谱
    gaussian_corr = exp(-(pi*fd*tau).^2);
end

% 相干时间计算 (50%相关性)
coherence_time_50 = 0.423 / fd;
% 相干时间 (90%相关性)
coherence_time_90 = 0.016 / fd;

%% 频率相关性分析
df = 0:10:2000; % 频率间隔 (Hz)

if use_parallel
    % 使用并行计算计算频率相关性
    fprintf('使用并行计算计算频率相关性...\n');
    parfor i = 1:length(df)
        % 指数型功率延迟分布
        freq_corr_exp(i) = 1/(1 + (2*pi*delay_spread*df(i))^2);
        % 均匀型功率延迟分布
        freq_corr_uniform(i) = sinc(delay_spread*df(i));
        % 高斯型功率延迟分布
        freq_corr_gaussian(i) = exp(-(pi*delay_spread*df(i))^2);
    end
else
    % 指数型功率延迟分布
    freq_corr_exp = 1./(1 + (2*pi*delay_spread*df).^2);
    % 均匀型功率延迟分布
    freq_corr_uniform = sinc(delay_spread*df);
    % 高斯型功率延迟分布
    freq_corr_gaussian = exp(-(pi*delay_spread*df).^2);
end

% 相干带宽计算 (90%相关性)
coherence_bandwidth_90 = 1 / (50*delay_spread);
% 相干带宽 (50%相关性)
coherence_bandwidth_50 = 1 / (5*delay_spread);

%% 封装结果
correlation_results = struct();
correlation_results.time_lags = tau;
correlation_results.jakes_corr = jakes_corr;
correlation_results.uniform_corr = uniform_corr;
correlation_results.gaussian_corr = gaussian_corr;
correlation_results.frequency_offsets = df;
correlation_results.freq_corr_exp = freq_corr_exp;
correlation_results.freq_corr_uniform = freq_corr_uniform;
correlation_results.freq_corr_gaussian = freq_corr_gaussian;

coherence_params = struct();
coherence_params.coherence_time_50 = coherence_time_50;
coherence_params.coherence_time_90 = coherence_time_90;
coherence_params.coherence_bandwidth_50 = coherence_bandwidth_50;
coherence_params.coherence_bandwidth_90 = coherence_bandwidth_90;
coherence_params.max_doppler = fd;
coherence_params.delay_spread = delay_spread;

%% 可视化结果
fig_correlation = figure('Name', '时频相关性分析', 'Position', [200, 200, 1200, 800]);

% 子图1: 时间相关性
subplot(2, 3, 1);
plot(tau*1000, jakes_corr, 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(tau*1000, uniform_corr, 'LineWidth', 2, 'Color', colors(2, :));
plot(tau*1000, gaussian_corr, 'LineWidth', 2, 'Color', colors(3, :));
plot([coherence_time_50*1000, coherence_time_50*1000], [-0.2, 0.5], '--', 'Color', colors(4, :));
plot([coherence_time_90*1000, coherence_time_90*1000], [-0.2, 0.9], '--', 'Color', colors(5, :));
grid on;
xlabel('时间间隔 (ms)');
ylabel('归一化自相关');
title('时间相关性');
legend({'Jakes模型', '均匀谱', '高斯谱', 'T_c (50%)', 'T_c (90%)'}, 'Location', 'NorthEast');

% 子图2: 频率相关性
subplot(2, 3, 2);
plot(df, freq_corr_exp, 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(df, freq_corr_uniform, 'LineWidth', 2, 'Color', colors(2, :));
plot(df, freq_corr_gaussian, 'LineWidth', 2, 'Color', colors(3, :));
plot([coherence_bandwidth_50, coherence_bandwidth_50], [0, 0.5], '--', 'Color', colors(4, :));
plot([coherence_bandwidth_90, coherence_bandwidth_90], [0, 0.9], '--', 'Color', colors(5, :));
grid on;
xlabel('频率间隔 (Hz)');
ylabel('归一化频域相关');
title('频率相关性');
legend({'指数型', '均匀型', '高斯型', 'B_c (50%)', 'B_c (90%)'}, 'Location', 'NorthEast');

% 子图3: 多普勒功率谱密度
subplot(2, 3, 3);
f_doppler = -fd:1:fd;
% Jakes谱
jakes_psd = 1./(pi*fd*sqrt(1 - (f_doppler/fd).^2));
jakes_psd(abs(f_doppler) >= fd) = 0;
% 均匀谱
uniform_psd = ones(size(f_doppler)) / (2*fd);
% 高斯谱
gaussian_psd = exp(-f_doppler.^2 / (2*(fd/2)^2));
gaussian_psd = gaussian_psd / sum(gaussian_psd) * (2*fd);

plot(f_doppler, jakes_psd, 'LineWidth', 2, 'Color', colors(1, :));
hold on;
plot(f_doppler, uniform_psd, 'LineWidth', 2, 'Color', colors(2, :));
plot(f_doppler, gaussian_psd, 'LineWidth', 2, 'Color', colors(3, :));
grid on;
xlabel('多普勒频率 (Hz)');
ylabel('功率谱密度');
title('多普勒功率谱密度');
legend({'Jakes谱', '均匀谱', '高斯谱'}, 'Location', 'NorthEast');

% 子图4: 相干时间与速度关系
subplot(2, 3, 4);
speeds = [1, 3, 10, 30, 100]; % 不同速度 (m/s)

if use_parallel
    % 并行计算相干时间
    fprintf('使用并行计算分析相干时间与速度关系...\n');
    parfor i = 1:length(speeds)
        coherence_times(i) = 0.423 / (speeds(i) * carrier_freq / c);
    end
else
    coherence_times = 0.423 ./ (speeds * carrier_freq / c);
end

semilogy(speeds, coherence_times*1000, 's-', 'LineWidth', 2, 'Color', colors(1, :));
grid on;
xlabel('移动速度 (m/s)');
ylabel('相干时间 (ms)');
title('相干时间与移动速度关系');

% 子图5: 相干带宽与时延扩展关系
subplot(2, 3, 5);
delay_spreads = [10e-9, 50e-9, 100e-9, 500e-9, 1e-6, 5e-6]; % 不同时延扩展

if use_parallel
    % 并行计算相干带宽
    fprintf('使用并行计算分析相干带宽与时延扩展关系...\n');
    parfor i = 1:length(delay_spreads)
        coherence_bws(i) = 1 / (50*delay_spreads(i));
    end
else
    coherence_bws = 1 ./ (50*delay_spreads);
end

loglog(delay_spreads*1e6, coherence_bws/1e6, 's-', 'LineWidth', 2, 'Color', colors(2, :));
grid on;
xlabel('时延扩展 (μs)');
ylabel('相干带宽 (MHz)');
title('相干带宽与时延扩展关系');

% 子图6: 时频相干区域
subplot(2, 3, 6);
% 创建相干区域图
time_range = [0.1, 1, 10, 100]; % 相干时间 (ms)
bw_range = [0.1, 1, 10, 100]; % 相干带宽 (MHz)

for i = 1:length(time_range)
    for j = 1:length(bw_range)
        % 计算对应的信道类型
        if time_range(i) > 10 && bw_range(j) < 1
            channel_type = sprintf('频率选择性\n慢衰落');
            color_idx = 1;
        elseif time_range(i) < 1 && bw_range(j) < 1
            channel_type = sprintf('频率选择性\n快衰落');
            color_idx = 2;
        elseif time_range(i) > 10 && bw_range(j) > 10
            channel_type = sprintf('平坦衰落\n慢衰落');
            color_idx = 3;
        else
            channel_type = sprintf('平坦衰落\n快衰落');
            color_idx = 4;
        end
        
        rectangle('Position', [log10(bw_range(j))-0.1, log10(time_range(i))-0.1, 0.2, 0.2], ...
                 'FaceColor', colors(color_idx, :), 'EdgeColor', 'k', 'LineWidth', 1);
        text(log10(bw_range(j)), log10(time_range(i)), channel_type, ...
             'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', ...
             'FontSize', 8, 'Color', 'w', 'FontWeight', 'bold');
    end
end
xlabel('相干带宽 (MHz)');
ylabel('相干时间 (ms)');
title('时频相干区域分类');
set(gca, 'XTick', log10([0.1, 1, 10, 100]), 'XTickLabel', {'0.1', '1', '10', '100'});
set(gca, 'YTick', log10([0.1, 1, 10, 100]), 'YTickLabel', {'0.1', '1', '10', '100'});

fprintf('时频相关性分析完成！\n');
fprintf('相干时间 (50%%): %.2f ms\n', coherence_time_50*1000);
fprintf('相干时间 (90%%): %.2f ms\n', coherence_time_90*1000);
fprintf('相干带宽 (50%%): %.2f MHz\n', coherence_bandwidth_50/1e6);
fprintf('相干带宽 (90%%): %.2f MHz\n', coherence_bandwidth_90/1e6);

end